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ABSTRACT 


Ray-theoretic  modeling  requires  accurate  amplitude  as  well  as  phase  both  for 
forward  modeling  and  Kirchhoff  inversion,  among  other  applications.  There  are 
no  analytical  solutions  to  the  ray  equations  in  realistic  earth  models,  thus,  we 
must  use  numerical  solutions  to  solve  problems  of  interest.  For  three  dimensional 
applications,  it  is  a  challenge  to  develop  numerical  modeling  codes  that  require 
reasonable  cpu  time  while  achieving  sufficient  amplitude  accuracy  to  be  useful 
in  applications. 

For  the  case  of  linear  sloth  (slowness  squared  or  inverse  wavespeed  squared), 
analytical  solutions  of  the  ray  equations  do  exist,  leading  to  a  combined  nu¬ 
merical  analytical  technique.  In  this  method,  the  physical  model  is  decomposed 
into  tetrahedral  blocks  of  sufficiently  small  size  to  allow  for  the  linear  sloth 
approximation  to  be  valid  in  each.  Analytical  solutions  in  each  tetrahedron 
are  then  pieced  together  to  provide  global  solutions.  The  ray  tracing  with  this 
method  is  relatively  fast.  However,  the  wavespeed  model  generated  by  this  tech¬ 
nique  is  not  sufficiently  smooth  to  produce  accurate  amplitudes,  numerically. 
Recent  attempts  to  further  smooth  the  physical  model  defeat  the  advantage  of 
speed  of  the  algorithm  because  the  smoothness  conditions  across  the  faces  of 
the  tetrahedra  generate  a  coupled  system  of  equations  of  a  size  proportional  to 
the  number  of  tetrahedra  in  the  global  physical  model.  This  is  not  practical  in 
3D,  Thus,  we  conclude  that  a  standard  smoothed  physical  model  on  a  Carte¬ 
sian  grid  is  likely  to  lead  to  a  computer  code  of  competitive  cpu  speed,  when 
amplitude  accuracy — dynamics — ^is  of  as  much  concern  as  traveltime  accuracy — 
kinematics. 

In  either  case,  we  use  a  wavefront  construction  technique,  in  which  the  size  of 
triangular  plates  connecting  three  nearby  rays  on  the  isochron  (surface  of  con¬ 
stant  traveltime)  is  used  as  an  indicator  of  adequate  density  of  rays.  When  the 
criteria  for  density  of  rays  are  violated,  data  at  new  points  on  the  wavefront 
are  interpolated  into  the  family  of  rays  and  the  wavefront  construction  contin¬ 
ues.  In  this  manner,  the  method  does  not  require  excessive  density  of  rays  at 
small  traveltimes  in  order  to  maintain  adequate  density  of  rays  at  larger  trav- 
eltimes.  The  technique  allows  for  multi-pathing  (caustics)  and  for  amplitude 
propagation  along  each  of  the  branches  of  the  wavefront. 

Applications  of  the  modeling  technique  are  shown. 

Key  words:  amplitude,  dynamic  ray  tracing,  analytic  ray  tracing,  wavefront 
construction 


Introduction 

In  this  report,  we  address  the  problem  of  accurate 
and  efficient  determination  of  multi-valued  3D  maps 
for  amplitude  as  well  as  traveltimes  or  ciny  other  ray- 
related  variables  throughout  the  target  zone  from  any 


shot  and  receiver  position.  The  current  interest  in  3D 
seismic  imaging  has  considerably  increased  the  impor¬ 
tance  of  ray  tracing  methods  in  wave  field  computa¬ 
tions.  Among  seismic  modeling  methods,  ray  tracing 
methods  provide  a  reasonable  compromise  between  ac- 
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curacy  and  computational  efficiency.  For  the  computa¬ 
tion  of  traveltime,  various  methods  have  been  described. 
Among  those,  the  finite-differencing  (FD)  method,  i.e., 
FD-solvers  of  the  eikonal  equation,  has  recently  become 
a  popular  method  for  calculation  of  “first  arrival”  travel- 
times  (Vidale,  1988).  However,  this  method  suffers  from 
the  disadvantages  that  it  is  restricted  to  the  compu¬ 
tation  of  first  arrivals  only  and  it  produces  unreliable 
amplitudes.  Both  are  severe  disadvantages  for  Kichhoff- 
type  algorithms,  such  as  the  Bleistein/Cohen  inversion 
(Bleistein  et  a/.,  1987;  Bleistein  et  al.^  1996),  where  the 
calculation  of  amplitudes  is  necessary  to  determine  the 
weighting  factor  in  inhomogeneous  media.  Furthermore, 
in  complex  media,  such  as  near  salt  domes  and  in  sub-salt 
regions,  later  arrival  traveltimes  should  be  considered  to 
obtain  better  image  quality.  Amplitudes  can  be  used, 
among  other  things,  to  find  most  energetic  arrivals. 

Simultaneous  computation  of  traveltimes  and  ampli¬ 
tudes  is  possible  by  dynamic  ray  tracing  (DRT).  It  pro¬ 
vides  accurate  multiple  arrivals,  amplitude  and  phase. 
Estimation  of  these  ray  data  can  be  carried  out  either 
by  numerical  solution  of  ray  tracing  equations  in  general 
smooth  grid-based  models  or  by  piecewise  analytic  solu¬ 
tions  for  certain  simple  velocity  functions  in  tetrahedral 
models.  Among  the  choices  of  ray  tracing  procedures, 
the  simplest  and  fastest  solution  of  the  ray  tracing  sys¬ 
tem  is  usually  based  on  its  analytic  solution,  wherever 
the  complexity  of  the  model  allows  one.  This  is  usually 
referred  to  as  analytic  ray  tracing  or  cell  ray  tracing. 
Generally,  the  whole  medium  is  divided  into  suitable  cells 
(usually  tetrahedra  in  3D),  in  which  the  velocity  can  be 
approximated  by  simple  functions  that  permit  analytic 
ray  solutions.  The  ray  in  the  whole  model  is  then  ob¬ 
tained  as  a  chain  of  analytically  computed  segments.  The 
analytic  ray  tracing  is  usually  performed  for  models  in 
which  either  the  velocity,  t;(xi),  or  l/t;(xi),  or  l/v^{xi), 
is  a  linear  function  of  Cartesian  coordinates.  The  sim¬ 
plest  analytic  solution  for  inhomogeneous  medium  is  the 
one  for  constant  gradient  of  squared  slowness,  also  re¬ 
ferred  as  linear  sloth  media  (Cerveny,  1987;  Meng  & 
Bleistein,  1997).  However,  this  assumption  leads  to  tetra¬ 
hedral  cells  with  artificial  second-order  discontinuities  at 
their  interfaces.  As  a  consequence,  this  approach  pro¬ 
duces  unreliable  amplitude  coefficients  across  the  inter¬ 
nal  boundaries.  Kdrnig  (1995)  proposed  a  method  us¬ 
ing  quadratic  sloth.  In  this  approach,  the  squared  slow¬ 
ness  and  its  gradient  with  respect  to  spatial  variables  are 
continuous  across  each  cell  boundary.  The  analytic  solu¬ 
tions  for  such  a  velocity  function  are  determined  by  using 
Laplace  transform.  Also,  the  computation  of  amplitude 
can  be  largely  simplified  by  calculating  the  ray  Jacobian 
directly  from  the  analytical  ray  equations.  However,  the 


problem  of  determining  the  cell  constants  in  quadratic 
sloth  is  rather  difficult.  The  model  design  leads  to  a  huge 
matrix  inversion  problem,  and  is  impractical  in  3D.  Only 
2D  implementation  of  the  traveltime  computation  was 
carried  out  by  Kdrnig  (1995).  We  have  concluded  that 
the  analytic  approach  in  tetrahedral  cells  does  not  likely 
offer  efficient  algorithms  in  dynamic  applications. 

Traditionally,  numerical  DRT  is  performed  by  shoot¬ 
ing  a  fan  of  rays  from  the  source  and  extrapolating  trav¬ 
eltimes  and  amplitude  away  from  the  rays  into  their 
nearby  regions  (Cerveny,  1987;  Virieux  Farra,  1991; 
Sun  &c  Biondi,  1995).  The  main  disadvantage  of  the  con¬ 
ventional  shooting  method  is  the  lack  of  control  of  ray 
density  in  the  search  fan.  Therefore,  it  is  hard  to  reach 
a  favorable  compromise  between  efficiency  and  reliabil¬ 
ity,  especially  in  complex  3D  models.  It  also  produces 
shadow  zones  in  areas  of  large  velocity  contrasts.  The 
wavefront  construction  traveltime  computation  method 
(Vinje  et  al.j  1996)  offers  a  solution  to  this  problem  by 
dynamically  adding  rays  as  needed.  In  this  method,  rays 
are  maintained  by  a  triangular  network,  and  are  traced 
stepwise  in  traveltime  through  the  model.  The  wave- 
fronts  are  then  obtained  automatically  as  a  by-product 
of  the  ray  tracing.  In  this  report,  the  idea  of  wavefront 
construction  is  applied  to  3D  complex  models  for  es¬ 
timation  of  both  traveltime  and  amplitude  coefficients. 
The  dynamic  interpolation  of  new  rays  assures  that  the 
wavefront  is  equipped  with  sufficient  ray  density  at  each 
computational  step.  Linear  interpolation  of  traveltime 
with  respect  to  the  simulated  wavefronts  and  linear  in¬ 
terpolation  of  amplitude  in  terms  of  tube  cross  sectional 
area  are  performed  at  grid  points  that  fall  into  the  sub¬ 
volume  formed  by  every  two  successive  wavefronts.  A 
grid  point  can  be  passed  by  different  sequences  of  wave- 
fronts  and,  thereby,  multi-valued  arrivals  can  be  detected 
and  recorded.  In  this  manner,  all  the  grid  points  in 
the  model  are  equipped  with  accurate-— perhaps  multi¬ 
valued — traveltimes  and  amplitudes. 

In  the  following  sections,  we  first  discuss  the  possi¬ 
bility  of  applying  analytic  solutions  in  tetrahedral  models 
for  amplitude  estimation.  Thereafter,  we  address  some 
important  issues  in  numerical  DRT  such  as  interpolation 
of  new  rays  and  estimation  of  parameters  at  grid  points. 
We  also  propose  a  smooth  gridded  model  representation 
for  the  purpose  of  computational  efficiency.  Finally,  we 
show  results  of  applying  this  method  to  different  velocity 
models. 

Dynamic  ray  tracing 

This  section  is  a  brief  review  of  the  dynamic  ray 
tracing  theory,  based  on  Cerveny  (1987),  (1995).  We  be¬ 
gin  by  introducing  two  coordinate  systems  involved  in 
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Figure  1.  The  ray-centered  coordinate  system,  ez  is  tangent 
to  the  ray  path  n.  el  and  el  form  a  plane  perpendicular  to  el, 
7i  and  72  are  the  ray  parameters  that  specify  the  ray,  usually 
they  are  either  the  take-off  angles(shown  in  this  figure)  or  the 
slowness  vector  components  at  the  source. 


DRT  —  the  ray-centered  coordinate  system  and  the  ray 
parameter  coordinate  system. 

The  ray-centered  coordinate  system,  denoted  by 
(gi,gf2,  ga),  is  a  curvilinear  orthogonal  coordinate  system 
associated  with  any  selected  ray  0  (see  Fig  1).  One  coor¬ 
dinate,  say  53,  corresponds  to  any  monotonic  parameter 
along  the  ray,  such  as  the  arc  length  s,  the  traveltime 
T  or  the  parameter  cr,  with  dr  =  dcr/v^.  Here,  we  take 
53  =  r,  the  traveltime  of  ray  Q.  away  from  the  source. 
Thus,  the  traveltime,  itself,  is  one  of  the  coordinate  axes 
in  the  ray-centered  coordinate  system.  Coordinates  51 
and  52  form  a  2D  Cartesian  coordinate  system  in  the 
plane  S  perpendicular  to  n  at  53  =  r,  with  the  origin  at 
f2.  The  vector  basis  of  the  ray-centered  coordinate  sys¬ 
tem  connected  with  0  is  formed  at  any  arbitrary  point 
^3  =  r  of  ray  O  by  a  right-handed  triplet  of  unit  vectors 
el(r),e2(r),e3(r),  as  shown  in  Fig  1,  Unit  vectors  ei(r) 
can  also  be  viewed  as  polarization  vectors  for  isotropic 
media.  The  unit  vector  el  determines  the  direction  of  the 
displacement  vector  of  P  waves,  which  is  always  linearly 
polarized.  Especially  important  are  unit  vectors  el, el, 
since  they  determine  the  polarization  of  S  waves,  when 
we  are  dealing  with  vector  solutions  of  the  elastic  wave 
equation. 

The  ray  parameter  coordinates,  (71,72,73),  are  de¬ 
fined  as  following:  71  and  72  are  the  ray  parameters  that 
specify  the  ray,  usually  they  are  either  the  take-off  an¬ 
gles  or  the  slowness  vector  components  at  the  source; 
73  is  any  monotonic  parameter  along  the  ray,  s,  r  or 
(7.  The  Jacobian  J  of  transformation  from  ray  coordi¬ 
nates,  (71,72,73),  to  the  general  Cartesian  coordinates, 
(xi,  X2,  X3),  is  an  important  factor  in  computation  of  the 
ray  amplitude  (Bleistein,  1984).  The  ray  amplitude  has 
the  following  form. 


A(x,xo)  = 


const 

W\' 


(1) 


with. 


J  =  det 


d{xi^X2i  X3) 

9(71,72,73) 


(2) 


Here,  A(x,xo)  is  the  amplitude  at  x  =  (xi,X2,®3)  cor¬ 
responding  to  the  source  xo  =  (xio,X2o,X3o);  J  is  the 
Jacobian;  and  the  constant  is  determined  by  the  choice 
of  (71,72,73)- 

The  rays  are  defined  as  the  characteristics  of  the 
eikonal  equation.  That  is,  by  transforming  the  eikonal 
equation  to  the  following  six  ray  equation  by  using  the 
method  of  characteristic,  (Bleistein,  1984) 


dxi 

dr 


v^Pi, 


dpi  _  1  dv 

dr  V  dxi  ’ 


1  =  1,2,3. 


(3) 


Here,  Xi  (r)  denotes  the  coordinates  of  position  along  the 
ray,  pi  (r)  denotes  the  components  of  the  slowness  vector, 
T  denotes  the  traveltime  along  the  ray,  and  v{xi)  denotes 
the  velocity.  System  (3)  is  often  referred  as  the  kinematic 
ray  tracing  (KRT)  system. 

Differentiating  the  KRT  system  (3)  with  respect  to 
the  ray  coordinates  7,  and  applying  Taylor  approxima¬ 
tion  up  to  second-order  in  qi  will  generate  the  dynamic 
ray  tracing  system.  The  DRT  system  can  be  expressed  in 
many  forms  and  in  various  coordinate  systems.  The  sim¬ 
plest  form  of  the  DRT  system  is  obtained  in  ray-centered 
coordinates  connected  with  the  ray  O: 


dr 


=  u^P, 


dr 


-“VQ, 

V 


where  Q,  P  and  V  are  2  x  2  matrices  defined  as 


Qij 

dqi  .  . 

=  hJ 

dyj 

=  1,2, 

Pij 

dpi  .  . 

=  h3 

d'yj 

=  1,2, 

V’’ 

d^V{qi,q2,S) 

dqidqj 

91=92=0 

=  Hkiv^kiHij, 

i,i  =  l,2, 

(4) 


Hki  =  ik'Ci.  (5) 

Here,  ik  are  the  basis  vectors  in  general  Cartesian  coor¬ 
dinates  (xi,X2,X3)  and  H  is  the  transform  matrix  from 
(91j92,53)  to  (xi,X2,a:3),  its  element  Hu  represents  the 
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/cth  Cartesian  component  of  basis  vector  e\.  Elements 
Vii  of  the  matrix  V  are  the  second  derivatives  of  veloc¬ 
ity  V  with  respect  to  qi,  and  equivalent  to  the  second 
derivatives  with  respect  to  Xi  under  transformation  H. 
Therefore,  dynamic  ray  tracing  (4)  requires  continuity  of 
the  velocity  field  up  to  second  derivatives.*  Q  is  a  trans¬ 
formation  matrix  from  the  ray  parameters  71 , 72  to  the 
ray-centered  coordinate  qi,q2-  Its  determinant  measures 
the  ray  Jacobian  (2),  and  is  also  called  the  geometrical 
spreading  factor.  P  is  a  transformation  matrix  from  the 
ray  parameters  71 , 72  to  the  slowness  vector  component 
in  the  ray-centered  coordinate  system.  The  quantity  P 
may  also  be  defined  by 


P  = 


KQ 


(6) 


where  K  is  the  wavefront  curvature.  The  KRT  system 
(3)  computes  the  first  derivatives  of  the  traveltime  field, 
while  the  DRT  system  (4)  relates  the  second  derivatives 
by  the  relationship, 


M(r)  =  PQ"‘. 


(7) 


where  M  is  2  x  2  matrix  M  of  second  derivatives  of  the 
traveltime  field  with  respect  to  the  ray-centered  coordi- 
nate  qi,q2,  Mi,j  =  d’^r/dqidqj. 

At  caustic  points,  the  ray  Jacobian  (2)  or  the  de¬ 
terminant  of  matrix  Q  vanishes.  In  3jD  structure,  there 
are  two  kinds  of  caustic  points  depending  on  the  range 
of  the  matrix.  The  ray  tube  may  shrink  to  a  caustic  sur¬ 
face  (envelope  of  rays)  which  is  perpendicular  to  the  di¬ 
rection  of  propagation  (a  caustic  point  of  the  first  or¬ 
der);  or  the  ray  tube  may  shrink  to  a  point  (a  caustic 
point  of  the  second  order).  In  passing  through  the  caus¬ 
tic  point  of  the  first  order,  the  ray  Jacobian  J  changes 
sign  and  the  argument  of  takes  on  the  phase  term 
Similarly,  in  passing  through  the  caustic  point  of 
the  second  order,  the  phase  term  is  ±7r.  The  phase  shift 
due  to  caustics  is  cumulative.  If  we  pass  through  several 
caustic  points  along  the  ray,  the  total  phase  shift  is  the 
sum  of  the  individual  phase  shifts,  this  is  often  referred 
as  the  KM  AH  amplitude. 


Analytic  ray  tracing 

It  is  known  that  the  realistic  velocity  field  of  interest 
is  often  rather  complicated  and  can  hardly  allow  a  gen¬ 
eral  analytic  solution  of  the  ray  tracing  system.  However, 
analytic  ray  tracing  plays  an  important  role  in  wave  field 
computation.  This  is  due,  in  part,  because  the  analytic 
solutions  are  valuable  in  the  cell  approach,  in  which  the 


whole  model  is  subdivided  into  a  set  of  tetrahedral  cells 
with  simple  velocity  functions  within  cells.  The  models 
allowing  analytic  solutions  are  usually  those  with  either 
the  velocity,  or  slowness,  or  squared  slowness  being  a  lin¬ 
ear  function  of  Cartesian  coordinates.  As  we  have  men¬ 
tioned  before,  this  constraint  does  not  provide  enough 
smoothness  for  amplitude  estimation.  In  this  section,  we 
discuss  the  analytic  solutions  for  quadratic  sloth  media. 

The  quadratic  sloth  distribution,  denoted  by  q,  is 
defined  as  a  quadratic  function  in  space, 

g(a;i,X2,X3)  =  -^7 - r 

v^{xi,X2,xz) 


=  A’{-2BiXi  +  CijXiXj.  (8) 

The  analytical  solutions  for  this  distribution  were  found 
by  Kornig  (1995)  using  the  Laplace  transform.  In  the 
Laplace  domain,  the  ray  coordinates,  Xi{s),  are  ratios 
of  polynomials  of  sixth  and  seventh  order,  respectively, 
in  s;  this  is  the  Laplace  variable  corresponding  to  cr,  the 
ray  tracing  integral  variable  with  da  =  v^dr.  The  ex¬ 
pressions  for  the  ray  trajectories,  Xi{a),  can  be  obtained 
explicitly  by  inverse  Laplace  transform  of  A'i(s)  using 
partial  fraction  expansions.  Depending  on  the  distribu¬ 
tion  of  eigenvalues  of  Cij,  the  solutions  of  Xi{a)  are  gen¬ 
eralized  into  seven  different  forms.  For  each  case,  the  ray 
trajectories  are  in  the  general  form 


Xi{a^  —  "^ikfk  (<^)7  ^  —  Ij2,3,  h  —  1,2,..., 7.  (9) 


Here,  the  wik^s  are  the  weighting  factors,  which  are  func¬ 
tions  of  Xio,pio,Bj  and  Cij  with  xio  and  pio  being  the 
initial  position  and  slowness  components;  fk{cr)  are  the 
basis  functions  corresponding  to  the  inverse  transform  of 
the  partial  fraction  expansions.  One  of  them  is  unity,  the 
others  are  either  low-order  polynomieils  in  a,  or  trigono¬ 
metric,  or  hyperbolic  function. 

Now  we  propose  an  alternative  to  Kornig’s  (1995) 
approach  to  calculate  the  amplitudes  along  rays.  Notice 
that  the  /a:(o‘)’s  in  (9)  are  functions  depending  only  on  <t, 
and  Wik^s  can  be  expressed  in  terms  of  XicPio,  and  con¬ 
stants  Bj  and  Cij .  Therefore,  if  we  choose  71  =  pio,  72  = 
P20  and  73  =  o-  in  (2),  the  ray  Jacobian  J  can  be  calcu¬ 
lated  analytically. 


J 


det 


feAw  ■ 

mk9k{cr)  . 


*  It  is  the  numerical  sampling  across  tetrahedral  interfaces  =  det{  (10) 

that  causes  amplitude  instability  when  the  second  derivative 

is  not  continuous.  with  summation  over  the  repeated  index,  k,  and 
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f  fki<r),  i  =  l,2 
\  9k{ff),  i  =  3. 


(11) 


The  coefficients  Jijk  can  be  expressed  for  all  the  seven 
cases  in  (9).  Notice  that  expression  (10)  has  taken  the 
place  of  the  DRT  system  (4),  and  it  is  less  computation¬ 
ally  costly  than  solving  the  eight  integral  equations  in 
small  time  steps. 

This  approach  of  making  quadratic  sloth  assump¬ 
tion  in  tetrahedral  models  has  eliminated  the  smoothing 
procedure  across  the  internal  interfcices.  However,  the 
constants  A,  Bi  and  Cij  in  (8)  are  usually  not  known  in 
advance  from  the  given  physical  model.  They  have  to  be 
determined  from  the  discrete  model.  The  assumption  of 
quadratic  sloth  is  equivalent  to  the  continuity  of  both 
q  and  the  gradient  of  q  across  the  internal  cell  faces. 
Therefore,  the  10  constants  A,  Bi  and  Cij  in  one  tetra¬ 
hedron  cannot  be  determined  by  the  velocity  values  at 
its  four  apices  only,  but  also  depend  on  the  values  in 
the  neighboring  cells.  Such  a  model  design  problem  for 
all  the  tetrahedral  cells  leads  to  a  huge  inverse  problem. 
If  the  whole  model  is  divided  into  N  cells,  the  size  of 
the  coupled  system  of  equation  is  proportional  to  ION, 
making  the  computation  very  time  consuming.  Further¬ 
more,  this  inverse  system  is  not  always  solvable,  or  has 
solutions  in  a  least  squares  sense,  at  best.  This  is  im¬ 
practical  in  3D  and  considerably  limits  the  applicability 
of  this  approach. 

Prom  the  above  discussion,  we  see  that  although  the 
assumption  of  quadratic  sloth  in  tetrahedral  models  pro¬ 
vides  accurate  amplitudes,  this  assumption  leads  to  a  dif¬ 
ficult  and  inefficient  numerical  problem  for  determining 
the  cell  constants.  And  this  problem  exists  in  all  exten¬ 
sions  to  quadratic  physical  models,  not  particularly  in 
a  quadratic  sloth  model.  We  conclude  that  analytic  ray 
tracing  in  tetrahedral  models  is  not  likely  to  give  us  an 
efficient  module  for  dynamics,  although  it  has  its  appli¬ 
cations  in  traveltime  calculations. 


Wavefront  construction  on  smooth  gridded 
models 

Here  and  below,  we  will  focus  on  numerically  solving 

(3)  and  (4)  for  smooth  gridded  models.  We  will  apply 
the  technique  of  wavefront  construction  (WFC)  to  both 
kinematics  and  dynamics. 

In  the  wavefront  construction  method,  a  relatively 
sparse  number  of  rays  are  shot  initially.  They  differ 
from  each  other  by  the  two  take-off  angles,  and  are  ex¬ 
trapolated  into  the  zone  of  interest  by  solving  (3)  and 

(4)  numerically  with  appropriate  initial  conditions.  Re¬ 
quired  accuracy  of  traveltime  and  amplitude  can  be  ap¬ 
proached  by  various  standard  numerical  procedures,  such 
as  Runge-Kutta  or  predictor-corrector,  for  example.  At 


any  computational  step,  the  wavefront  is  obtained  as 
a  by-product  of  the  ray  tracing.  The  wavefront  is  rep¬ 
resented  by  triangular  plates  that  connect  every  three 
neighboring  ray  endpoints  on  the  WF.  The  nearby  rays 
in  3D  are  then  defined  and  organized  by  such  a  triangu¬ 
lar  mesh  consisting  of  the  internal  ordering  of  connect¬ 
ing  endpoints  in  each  triangle  and  adjacent  triangle (s)  to 
each  of  its  sides.  The  processes  of  checking,  interpolat¬ 
ing  of  new  rays  and  estimating  of  grid  point  parameters 
(described  below)  are  all  performed  within  such  a  trian¬ 
gular  network.  Rays  are  added  and  the  original  triangle 
is  subdivided  into  new  triangles  when  certain  criteria, 
restricting  the  size  of  the  triangular  plates,  are  violated. 
In  this  manner,  the  wavefront  always  has  sufficient  ray 
density  without  a  priori  estimation  of  the  number  of  rays 
needed,  or  by  imposing  an  excessive  ray  density  on  ini¬ 
tiation.  For  complex  3D  velocity  models,  the  wavefront 
surface  may  be  very  complicated,  folding  in  on  itself  at 
some  parts,  for  example;  however,  no  tears  or  holes  in 
the  interior  of  the  surface  are  allowed.  In  this  sense,  the 
wavefront  is  complete.  On  the  other  hand,  a  grid  point 
can  be  passed  by  different  sequences  of  wavefronts;  multi¬ 
valued  arrivals  can  then  be  estimated  and  distinguished 
by  their  initial  take-off  angles. 

The  most  attractive  advantage  of  the  WFC  method 
is  that  it  is  more  efficient  than  the  conventional  ray  trac¬ 
ing  method.  In  addition,  WFC  gives  better  ray  coverage, 
especially  in  areas  of  large  geometrical  spreading  where 
conventional  ray  tracing  may  give  no  arrivals.  Further¬ 
more,  compared  to  FD-solvers,  the  WFC  method  is  not 
restricted  to  the  calculation  of  only  first  arrivals.  Ampli¬ 
tude  and  other  ray  theoretical  quantities  are  also  avail¬ 
able.  Thus,  it  meets  the  requirements  for  accurate  mod¬ 
eling  of  amplitude  as  well  as  phase,  a  requirement  for 
inversion  as  opposed  to  migration. 

Ray  interpolation 

The  wavefront  construction  method  is  largely  de¬ 
pendent  on  the  procedure  of  interpolation  of  the  wave- 
front  at  each  step.  New  ray  endpoints  must  be  added 
along  the  simulated  wavefront  and  must  have  the  prop¬ 
agation  direction  that  the  ray  would  have  had  if  it  had 
been  shoot  from  the  source.  This  section  addresses  an 
algorithm  for  this  procedure. 

Rays  diverge  and  the  wavefronts  expand  through  the 
wave-field.  When  new  ray  endpoints  are  needed  to  keep 
a  certain  ray  density  on  the  wavefront,  the  whole  trian¬ 
gular  network  will  have  to  be  reorganized.  The  criterion 
for  this  interpolation  can  be  the  area  of  triangles,  which 
must  not  exceed  a  pre-specified  limit,  and/or  the  angle 
deviation  of  the  slowness  vectors  of  two  adjacent  rays, 
which  cannot  be  too  large.  New  rays  are  always  added 
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Figure  2.  Interpolation  of  new  rays  and  estimation  of  grid  parameters  are  performed  in  a  “ray  tube”.  (a)Interpolation  of  a 
single  new  ray.  The  two  ray  endpoints  and  their  propagation  directions  form  two  straight  lines  in  3D.  An  “approximate”  center 
can  be  defined  as  the  midpoint  of  the  line  segment  that  connects  the  two  straight  lines  at  their  points  of  shortest  distance.  This 
approximate  center,  along  with  the  two  ray  ends  at  the  old  WF,  form  a  fan  and  a  circular  curve  connecting  the  two  ray  end 
points.  The  new  ray  position  is  then  found  along  the  dividing  direction  from  the  approximate  center,  and  at  its  intersection 
with  the  circular  curve.  Then  the  interpolated  ray  is  traced  from  the  old  to  the  new  WF.  (b) Simple  ray  cell  with  an  interior 
grid  point.  Ray  data  are  estimated  with  respected  to  the  two  simulated  WFs. 


Figure  3.  3D  wave  field  of  a  linear  sloth  model  using  WF 
construction  method.  The  grey  part  is  the  shadow  zone. 


in  between  the  pairs  of  existing  rays  in  order  to  meet  the 
criteria  of  size  and  shape  of  triangular  plates  on  the  wave- 
front.  Fig  2-a  illustrates  the  interpolation  of  a  single  new 
ray.  The  two  ray  endpoints  and  their  propagation  direc¬ 
tions  form  two  straight  lines  in  3D.  An  “approximate” 
center  can  be  defined  as  the  midpoint  of  the  line  seg¬ 
ment  that  connects  the  two  straight  lines  at  their  points 
of  shortest  distance.  This  approximate  center,  along  with 
the  two  ray  ends  at  the  old  WF,  form  a  fan  and  a  circular 
curve  connecting  the  two  ray  end  points.  The  new  ray 
position  is  then  found  along  the  dividing  direction  from 
the  approximate  center,  and  at  its  intersection  with  the 


Figure  4.  The  rays  (white)  of  linear  sloth  model  have  a  para¬ 
bolic  shape.  The  grey  surface  is  the  caustic  surface  of  the  ray 
equations  for  this  model. 


circular  curve.  Other  parameters  along  the  new  ray  are 
interpolated  linearly. 

There  are  other  alternative  methods  of  interpolating 
new  rays,  such  as  the  parameterization  of  a  wavefront  by 
a  third-order  polynomial  (Vinje  et  al.^  1996).  However, 
they  require  special  treatment  in  the  vicinity  of  caustic 
points,  since  only  rays  belonging  to  the  same  phase  must 
be  used  to  determine  quantities  of  the  new  ray.  Moreover, 
we  do  not  use  the  curvature  of  the  wavefront  obtained 
from  DRT  for  the  interpolation  to  keep  the  problems  of 
KRT  and  DRT  separate.  We  have  found  this  method 
very  stable. 
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Grid  interpolation 

Another  interpolation  procedure  in  DRT  is  the  es¬ 
timation  of  ray  data  at  the  output  grid  points.  This  is 
our  final  goal  of  the  ray  tracing  algorithm.  We  perform 
the  grid  interpolation  within  ray  tubes,  which  are  prism¬ 
shaped  bodies  bounded  by  three  rays  and  the  triangles 
that  connect  them  on  the  two  WFs.  First  the  grid  points 
falling  into  (or  close  to)  each  cell  are  found.  Then,  the 
traveltime  can  be  estimated  at  the  each  grid  point  in 
a  similar  way  to  the  interpolation  procedure  for  new 
rays(see  Fig  2-b).  The  approximate  center  is  determined 
by  the  three  rays  with  ray  endpoints  on  any  of  the  two 
WFs.  The  distances  from  each  grid  point  to  the  approx¬ 
imate  center  and  to  the  simulated  wavefront  are  calcu¬ 
lated.  The  traveltime  at  the  grid  point  is  then  recorded 
as  to  +  d/u,  with  to  the  time  at  the  wavefront,  d  the  dis¬ 
tance  of  the  grid  point  away  from  the  wavefront,  and  v 
the  velocity  at  the  grid  point. 

The  above  procedure  is  not  suitable  for  interpolation 
of  amplitude,  because  the  isochrons  surface  is  usually  not 
the  iso-amplitude  surface.  Here,  we  apply  linear  inter¬ 
polation  for  amplitudes,  which  is  based  on  the  assump¬ 
tion  that  the  amplitudes  vary  only  slowly;  otherwise,  the 
validity  conditions  of  the  underlying  asymptotic  theory 
would  be  violated.  Since  the  ray  Jacobian — the  determi¬ 
nant  of  Q  in  (4) — is  proportional  to  the  cross  sectional 
area  of  the  ray  tube,  we  interpolate  the  ray  Jacobian 
linearly  with  respect  to  the  triangle  areas  on  the  two 
wavefronts. 

Another  parameter  that  requires  interpolation  at  ev¬ 
ery  grid  point  is  the  initial  shooting  direction,  i.e.,  the 
initial  take-off  angles  of  the  ray  that  would  reach  the 
grid  point  if  the  ray  actually  had  been  traced.  This  pa¬ 
rameter  is  stored  in  order  to  distinguish  between  arrivals 
because  two  arrivals  at  a  grid  point  cannot  have  almost 
equal  take-off  angles. 


Model  representation 

The  smoothness  of  the  velocity  model  representation  is 
critical  to  the  calculation  of  amplitudes.  The  integration 
of  the  DRT  system  (4)  requires  continuity  of  the  veloc¬ 
ity  field  up  to  the  second  derivatives.  Many  ray  tracing 
procedures  (Farra,  1990)  involve  a  type  of  spline  interpo¬ 
lation  for  the  evaluation  of  velocities  at  arbitrary  points. 
Spline  interpolation,  however,  is  a  time  consuming  pro¬ 
cedure.  Here,  we  define  the  velocity  model  on  a  fine  grid 
(about  three  or  four  grid  points  per  shortest  significant 
model  wavelength)  and  pre-calculate  its  first  and  second 
derivatives  at  all  grid  points  by  finite  differences  of  sec¬ 
ond  order.  For  the  evaluation  of  the  velocities  and  their 
derivatives  at  arbitrary  points  we  use  linear  interpola- 
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Figure  5.  Four  cross  sections  of  the  3D  SEG  saltdome  veloc¬ 
ity  model  (subsection)  at  i  =  1.5fcm,a:  =  S.5km,y  =  1.5A;m, 
and  y  =  3.5fcm. 
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Figure  6.  Triangulated  wavefronts  for  smoothed  salt  dome 
model  at  t  =  0.1,0.4,0.65,  and  0.9s.  The  source  is  located  at 
the  center  of  the  up  surface. 

tion.  For  the  smooth  models  defined  on  fine  grids,  the 
difference  between  this  linear  representation  and  a  spline 
representation  of  the  model  is  negligible. 

When  the  considered  model  contains  discontinuous 
velocities,  a  smoothing  procedure  must  be  applied  to 
guarantee  that  the  velocities  vary  smoothly.  For  the  sake 
of  computational  efficiency,  the  interface  conditions  axe 
eliminated  here  by  using  proper  smoothing  and  a  densely 
sampled  grid  model. 

Examples 

The  first  exzunple  provides  a  test  of  the  accuracy  of 
this  modeling  technique.  The  synthetic  model  we  choose 
for  this  test  is  the  one  with  constant  gradient  of  squared 
slowness,  i.e., 

'’■57-r  =  ^  +  BiXi  H-  B2X2  H-  Bzxz.  (12) 

u2(x) 

In  such  a  medium,  the  rays  have  a  parabolic  shape;  both 
traveltime  and  amplitude  field  can  be  expressed  exactly 
by  analytical  solutions  for  comparison  purposes.  Figure 
3  shows  the  3D  wave  field  with  the  gradient  constants 
being  (0, 0,  —0.2)s^/km^.  The  relative  difference  of  com¬ 
puted  traveltimes  and  analytic  ones  is  less  than  0.1% 
through  the  whole  interest  area,  while  the  differences  be¬ 
tween  computed  and  analytical  amplitudes  are  no  more 
than  1%.  The  grey  part  in  Fig  3  is  the  shadow  zone, 
where  no  rays  are  entered.  This  is  due  to  the  negative 
gradient  constant  of  Bz,  which  is  equivalent  to  the  in¬ 
creasing  velocity  with  depth.  In  Figure  4,  the  grey  sur¬ 
face  is  the  caustic  surface  of  the  ray  equations.  It  is  the 
envelope  of  the  parabolic  rays. 

The  aim  of  the  second  example  is  to  prove  that 
our  tracing  algorithm  can  operate  on  a  complex  velocity 
model.  This  example  is  performed  on  a  3D  SEG  salt  dome 


velocity  model.  This  synthetic  velocity  model  contains 
one  complex  salt  dome  with  high  velocity  in  the  dome 
and  low  and  slowly  varying  lateral  layers  outside  the  salt 
dome.  Due  to  the  RAM  capacity  of  the  computer,  we 
extracted  a  subsection  of  the  original  model,  which  has 
part  of  the  salt  in  the  middle.  The  strong  velocity  con¬ 
trast  at  the  salt  wall  has  violated  the  smoothness  re¬ 
quirement  of  ray  tracing,  therefore  we  first  smoothed 
this  reference  model.  Figure  5  shows  four  cross  sec¬ 
tions  of  the  smoothed  velocity  model.  The  grid  size  is 
40m  X  40m  x  40m.  Figure  6  shows  some  wavefronts  for 
this  model.  The  wavefronts  expand  and  become  more 
complex  for  the  later  traveltimes.  However,  by  using 
the  wavefront  construction  technique,  all  the  wavefronts 
have  sufficient  ray  density.  Figure  7  presents  cross  sec¬ 
tions  of  traveltime  maps  of  the  above  smoothed  velocity 
model.  The  isochrons  are  spherical-like  at  the  shsillow 
parts,  but  the  shape  changes  due  to  the  salt  in  the  mid¬ 
dle  depth. 

Conclusions 

We  have  demonstrated  that  numerically  solving 
the  ray  equations  on  a  smooth  gridded  model  pro¬ 
vides  forward  modeling  that  is  fast  enough  for  three- 
dimensional  computations.  Both  the  traveltimes  and  am¬ 
plitudes  proved  to  be  smooth  and  stable  in  our  examples. 
It  remains  to  check  the  numerical  accuracy  of  this  tech¬ 
nique,  as  compared  to  analytical  solutions  and  alterna¬ 
tive  numerical  methods.  However,  it  is  already  known 
that  the  tetrahedral-based  approach  produces  unaccept¬ 
able  amplitudes,  due  to  the  difficulty  in  efficiently  obtain¬ 
ing  accurate  amplitudes  across  the  internal  interfaces. 

The  WFC  procedure  based  on  proper  interpolation 
of  new  rays  makes  the  dynamic  ray  tracing  more  effi¬ 
cient  and  results  in  a  dense  and  consistent  ray  coverage 
throughout  the  model,  even  in  areas  of  large  geometriccd 
spreading.  When  accurate  amplitudes  are  required,  we 
believe  that  this  is  a  competitive  method  for  develop¬ 
ment  of  ray  theoretic  wavefields. 
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